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^ ' We study the linear stability of Plane Poiseuille flow of an elastoviscoplastic 

^N , fluid using a revised version of the model proposed by Putz and Burghelea 

(Rheol. Acta (2009)48:673-689). The evolution of the microstructure upon a 
P\J ■ gradual increase of the external forcing is governed by a structural variable 

f~^ . (the concentration of solid material elements) which decays smoothly from 

^D I unity to zero as the stresses are gradually increased beyond the yield point. 

Stability results are in close conformity with the ones of a pseudo-plastic fluid. 

Destabilizing effects are related to the presence of an intermediate transition 

zone where elastic solid elements coexist with fluid elements. This region 
^ ' brings an elastic contribution which does modify the stability of the flow. 

1. Introduction 

During the past several decades, yield stress fluids found an increasing number of prac- 
tical applications for several major industries (which include cosmetics, foods, oil field 
etc.) and they are encountered in the daily life in various forms such as hair gels, food 
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pastes, cement, mud and more. Such materials can be loosely and simplistically de- 
fined as materials that do not flow unless a minimal stress (referred as 'yield stress') is 
applied onto them. Understanding and controlling the hydrodynamic stability during 
flows of yield stress fluids is important for practical applications which involve flows of 
such materials through conduits. From a fundamental standpoint, yield stress materials 
continue triggering intensive debates and posing difficult challenges, of both theoretical 
and experimental nature. Undoubtedly, the best known and most cited debate is related 
to the very definition of such materials and the existence of a 'true yield stress', [1,2]. 
It was argued in Refs. [1,2] that the yield stress emerges as an artefact related to the 
inability of the rheometric equipment to properly identify a viscous flow regime in a range 
of negligibly small rates of deformation. Recent studies claim to have solved the "yield 
stress debate" by arguing that prior to yielding the viscosity of the material is infinite 
which demonstrates the existence of a true solid state, [4,23]. 

1.1. Yielding of a Carbopol® gel and its relevance to hydrodynamic 
studies 

The existence of a "true" yield stress and the nature of the transition from a solid to fluid 
behaviour may look at a first glance of little or no relevance to the hydrodynamic stability 
of the flow of a yield stress material. Indeed, whereas the yielding transition occurs at low 
values of the Reynolds number (typically Re < 1) a loss of the hydrodynamic stability is 
typically observed at significantly larger Re (typically Re > 1000). On the other hand, 
the base flows usually considered in the linear analysis of the hydrodynamic stability of 
channel flows of yield stress fluids are characterized by a significant stratification of the 
velocity gradients: large values near the channel boundaries which are consistent with a 
yielded flow region and vanishing values near the center-line, which are consistent with 
a plug region. A recent experimental investigation of the laminar-turbulent transition 
in the pipe flow of a yield stress fluids demonstrates that the transition to turbulence 
occurs when the Reynolds stresses balance the yield stress of the fluid, that is when the 
solid plug is broken, [14]. These findings corroborate well with the idea that, in fact, the 
nature of the solid-fluid transition and the yielding scenario may play an important role 
in the stability problem. A rather new and certainly unexpected insight into the yielding 
transition has been brought by an experimental study on sedimentation of spherical par- 
ticles in a physical gel (Carbopol® 940), [35]. Quite unexpectedly, the fore-aft symmetry 
of the flow patterns around the falling object was broken (in spite of the smallness of 
the Reynolds number during these experiments) and even more surprisingly a negative 
wake has been observed (see Fig. 8 in Ref. [35]). A similar non symmetric flow pattern 
has been observed experimentally during flows of Carbopol® gel past a cylinder, [41]. A 
first message of the sedimentation study presented in [35] was that the solid-fluid tran- 
sition in the Carbopol® gel might not be reversible upon increasing/decreasing stresses, 
unlike it has been traditionally assessed. A second important message was that the elas- 
tic effects responsible for the emergence of the experimentally observed negative wake 
might be consistent with an intermediate transition regime where solid material elements 
(unyielded) coexist with fluid (yielded) ones. A more systematic investigation of these 



features of the solid- fluid transition in a Carbopol® gel has been presented in Ref. [34]. 
Indeed, as suggested by the sedimentation experiments presented in [35], a hysteresis has 
been found in the solid-fluid transition of a Carbopol® gel upon increasing/decreasing 
stresses (see Fig. 3 in Ref. [34]). A clear signature of the elasticity involved in the tran- 
sition region in the form of a recoil effect (negative shear rates) on the decreasing stress 
branch of the flow curve has been observed as well. The hysteresis effect observed in 
the increasing/decreasing stress flow curves and reported in [34] has also been observed 
experimentally by Divoux et al., [9], and attributed to critical slowing down of the system 
dynamics close to yielding (and not to any memory effect). 

1.2. Linear stability analysis of a yield stress fluid: a brief survey 

The first study on linear stability analysis of a Bingham viscoplatic fluid is due to Frigaard 
et al. [10]. These authors concentrated on the asymptotic stability of two-dimensional 
traveling wave perturbations and found that the critical Reynolds number increases lin- 
early with Bingham number. They also showed that the plug region remains unaffected 
by the disturbance field. In this work the authors imposed even symmetry at the yield 
surface and the results are incomplete. Some years later Nouar et. al. studied again 
this problem using the exact boundary conditions at the yield surface and they found 
that the plane-Bingham-Poiseuille fiow is linearly stable. This being a consequence of 
the vanishing perturbation at the yield surface. They also conjectured that this result 
could be extended to Poiseuille flow of a viscoplastic fluid in a circular or annular pipe 
and to other rheological models such as the Herschel-Bulkley and Casson models. 

As mentioned before, recent studies have found that the solid-fluid transition in some 
yield stress materials, such as Carbopol®, is not reversible under increasing/decreasing 
stresses, thus the objective of this work is the linear stability analysis of an elastovis- 
coplastic fluid. In this work we consider a modified version of the structural model devel- 
oped by Putz and Burghelea in [34]. This models consists of a non- linear Maxwell- type 
viscoelastic constitutive equation and a kinematic equation that governs the behaviour 
of the microstructure. The relaxation time of the fluid depends on the structural variable 
which is equal to one if the fluid is unyielded and zero if the fluid is fully yielded. We 
should note that the transition from viscoelastic solid to fluid is continuous and smooth. 
For the non-linear viscosity we consider a regularized viscoplastic model, its behaviour 
is the one of a pseudo-plastic fluid, that is, if the second invariant of the stress tensor is 
below the yield stress the fluid presents a very high viscosity. 

Several authors have studied the linear stability problem of a flow with variations in 
viscosity across the channel [6,25,32]. Govindarajan et. al. [13] showed that by having 
a viscosity which is decreasing function of space it is sufficient to considerably delay 
the onset of two-dimensional instabilities. This result was latter confirmed by Nouar and 
Bottaro for the case of shear-thinning fluids [25] . Linear stability analysis for a regularized 
Bingham model were carried out by Frigaard and Nouar in [11]. They showed that the 
spectrum for the regularized problem had some physical spurious eigenvalues, not present 
in the Bingham problem. The stability of these eigenvalues depended on the size of the 
regularization parameter and unstable modes can be found even for moderate Reynolds 



numbers if this parameter is small enough. 

For the case of Maxwell liquids, Gorodstov and Leonov [12] showed the existence 
of a stable continuous spectrum and two stable eigenvalues for each value of the wave 
number in the streamwise direction, a. In the same work they predicted instabilities 
at finite Reynolds numbers. Later, this result was proven wrong by several authors, 
e.g. [16, 20, 36] . Through careful numerical simulation of the generalised eigenvalue 
problem no unstable modes were found, even for high Reynolds numbers. Denn and 
coworkers [16,33,38] also studied the stability of plane Poiseuille flow and found that at 
high Reynolds numbers UCM fluids are significantly less stable than Newtonian fluids. 
Sureshkumar & Beris [40] confirmed this and also showed that destabilisation is reduced 
for the Oldroyd-B and Chilcott-Rallison models. Renardy [37] has proven the linear 
stability of Couette flow of a UCM fluid, but there is no proof of linear stability of 
Couette flows for more general fluids of this type. However, it is generally believed that 
plane Couette flow is linearly stable. Indeed much of the literature is focused at the study 
of interesting features of the eigenspectra, rather than marginal stability, e.g. [19,42]. 

2. Elastoviscoplastic model with internal microstructure 

The elastovisplastic fluid is described by the following set of equations with unknowns 
(p,-u,T,$): 

f du \ 

p[^ + {u-V)u\ = -Vp + Hs'V -j + V -T (1) 

V-u = 0, (2) 

where fig denotes the solvent viscosity. The rate of strain tensor is 7 = 'y{u). The 
constitutive equation for the elastic stress tensor r is: 

, + ^i(^<|,^=,(^H)^, (3) 

where 

• = Vw Vw^ 

Dt 

is the Upper-Convected-Derivative and D ■ /Dt is the usual material derivative. Finally, 
the concentration of the solid state, $, satisfies the following kinematic equation: 

— + {u-V)^ = Rd{^, t{u)) + Rr{<^, t{u)), (4) 

with Rd being the rate of destruction of solid units and Rr is the rate of fluid recombi- 
nation of fluid elements into a gelled structure. 

Note that when <I> = 1 the model describes a Maxwell type viscoelastic fluid with a 
nonlinear viscosity nij) and if <I> = the model represents a Generalized Newtonian 
fluid. 




In both (3) and (4), 7 and r are the second invariant of the rate of strain and elastic 
stress tensors respectively. These are defined as follow: 

3 ^ ^/^ 

2.1. Choice of functions for nonlinear viscosity fi^j), destruction rate 

Rd{^,r{u)) and recombination rate Rr{^,T{u)). 

In this section we present our choices for the nonlinear viscosity ^(7) and the destruction 
and recombination rate functions. The choices are clearly non unique. Constitutive 
models for fluids with a yield stress trace back to the early 1900's . The most widely 
used are the ones due to Bingham [3], Herschel-Bulkley [15] and Casson [5]. Description 
of the rates of destruction and recombination is rather empirical and for this we follow 
closely [34]. 

2.1.1. Non-linear viscosity fi{'j) 

The fluids we are interested in present a yield stress. It is well known that an "ideal" 
viscoplastic fluid behaves as a Newtonian (or Generalized Newtonian) fluid if r is above 
the yield stress, if the contrary happens then the material behaves as a plastic solid 
and one says that the fluid is unyielded. In [34] it has been shown that the transition 
from an elastic solid behavior (unyielded) to a fluid behavior (yielded) is not direct but 
mediated by a viscoelastic like regime where solid and fluid behavior coexist. This effect, 
though previously unnoticed, is particularly pronounced in the case of "fast" yielding, 
i.e when the externally applied stress changes unsteadily and becomes vanishingly small 
in the limit of steady yielding, in other words, when the stresses vary infinitesimally 
slow. Although it has been argued that if one waits "long enough" the viscosity in the 
unyielded region actually tends to infinity, see [23], in an experimental (to investigate 
hydrodynamic stability for example) or industrial setting it is unlikely one will reach those 
time frames. Therefore for the rest of this study we will work with a regularized version of 
a viscoplastic constitutive model. Because our main interest is on linear stability analysis 
and this involves high Reynolds numbers it is better if we have a model with non-zero 
infinite shear viscosity, for this reason we choose the shear-thinning regularized Casson's 
model for the viscosity: 
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where m > 1 is the power law index and Ty is the yield stress and e <C 1 is the regu- 
larization parameter . We recover Casson's model when m = 2, Bingham model when 
m = 1 and the model reduces to a Newtonian viscosity when m = 1 and Ty = 0. In 
Figure la show an example of (5), note that if r < r^ the viscosity becomes very large 



but never infinity as long as e is fixed and different from zero, thus we will be working 
with a pseudo-plastic viscosity (shear-thinnig). 

2.1.2. Rdi<^,T{u)) and Rri^,T{u)) 

As discussed by Putz and Burghelea in [34] Rd and Rr have to satisfy the following 
assumptions: 

1. Rdi^, 't{u)) is proportional to the relative speed of neighboring solid units and the 
existing amount of solid, thus 

Rd{^,T{u)) = -g{T{u))^, (6) 

and 

giT{u)) = K, (l + tanh (^^^T^) ) ' ^^^ 

where g{T) is an increasing function of the second invariant of the stress tensor, 
Kd is the rate of destruction with units s~^ and w determines the "width" of the 
solid-fluid region and has units Pa. 

2. Rr{^, ^i^)) decreases with the relative speed of neighboring elements, begin prac- 
tically zero in a fast enough flow. Therefore, 

Rr{^,T{u)) = f{T{u)){l-'^), (8) 

and 

/(rH) = Kr (l -tanh (^^^T^)) ' ^^^ 

where /(r) is a decreasing function of the second invariant of the stress tensor. In Figure 
lb we show an example of (7) and (9). Note that when r is far below the yield stress, 
the material is fully "unyielded" which means it behaves as a viscoelastic solid. This 
because when r < Ty, nij) is of the order 0(l/e) and <l> = 1, thus by a simple dominant 
balance analysis equation (3) reduces to the one for a Kelvin- Voight viscoelastic solid. 
If the opposite happens (r far above the yield stress) $ = and the fluid behaves as 
a shear-thinning fluid. It is important to note that in a neighborhood around Ty the 
material is neither a viscoelastic solid nor a viscous fluid, both solid and fluid structures 
coexist and this region will play a crucial role int the stability of the flow. The existence 
of this region is the main difference with respect to other proposed models for this type 
of elastoviscoplastic fluids. The main difference of the approach used here and previous 
work, see [22], is that we do not use the steady states of (4) but solve it numerically and 
thus describe a non-steady yielding behavior which we believe might be of some practical 
relevance 
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Figure 1: (a) Viscosity function /i(7) with m = 3, ^oo = 0.0102, Ty = 6.5 and e = 0.01. 
(b) /(r) and ^(r) with i^d = i^,. = 0.3 and w = 0.1. 



2.2. Validation of the model against experimental results 

To validate the model, we use the rheological measurements presented in [34] conducted 
for various Carbopol® solutions with concentrations (by weight) ranging in between 
0.05% and 0.2%. The measurements were conducted on a Bohlin (now Malvern) C—VOR 
rotational rheometer equipped with a digitally controlled temperature (within 0.1° C ) 
bath [RTE — 111 from Neslab). The rheometer was operated in a controlled stress 
mode. To insure a good reproducibility of the measurements, prior to each experiment the 
sample has been pre-sheared at a constant stress (usually the largest stress applied during 
the test, so that yielding occurs in the entire volume of the material) for 300 s. After 
the pre-shear step, the samples were allowed to relax for 300 s prior to the rheological 
measurements. A first major concern during these experiments was the occurrence of the 
wall slip phenomenon at the contact with the measuring geometry, which is a well known 
and documented effect for Carbopol® gels, [31]. To prevent this, a serrated parallel 
plates geometry with an approximate roughness of 0.8 mm has been used. A second 
concern was related to the possible artefacts introduced by fluid evaporation during long 
experimental runs. To prevent this, a solvent trap has been placed around the free fluid 
meniscus. After each experimental run we have carefully checked (by visual inspection) 
that no significant changes in the shape of the meniscus occurred, and thus concluded 
that evaporation effects were either minimal or absent. The radius of the plates was 
R = 2 cm. and the distance between them was d = 1 mm (measured between the plates 
protuberances). Two types of rheological tests have been performed: controlled stress 
linear ramps and controlled stress oscillatory sweeps at a fixed frequency. For each sample 
we have conducted controlled stress experiments for 19 different values of the total ramp 
time. Each constant stress rheological experiment started with a increasing stress ramp 
and ended with a decreasing stress ramp within the same range of stresses and the 
same stress step. The data averaging time per stress value, Iq (referred in [34] as the 
characteristic time of forcing), has been varied between 0.2 s and 2 s. For each up-down 
stress ramps 1000 stress values have been explored ranging in between 0.1 Pa and 20 Pa. 
We emphasize here that a true steady state of deformation can only be inferred in the 
asymptotic limit to ^ oo, unlike in previous studies concerning Carbopol® gels where 
a steady state was a priori set by deliberately choosing large values of to (an accurate 
description of such procedure is presented, for example, in Ref. [7]). Alternatively, the 
time dependent response of the samples has been tested via stress controlled oscillatory 
experiments at several frequencies and in the same range of stresses. 

We now turn our attention to the comparison of the predictions of the model described 
in Sec. 2 with the experimental results presented in Ref. [34]. 

Bearing in mind that the experiments provide a scalar data set, the equations (l)-(4) 
reduce to: 



^ = -ff(T)<I> + /(r)(l-c|>) (10) 

^<&§+r = M7)7 (11) 



Because the rheometer operates in a controlled stress mode, the above equations de- 
couple. First we solve for (10) using the function odelSs provided with the MATLAB® 
distribution. The function odeISs performs well with stiff problems and its accuracy is 
recognized as medium to high. Equation (11) reduces to an algebraic equation for 7 
which we solve using the function fzero provided with the MATLAB® distribution as 
well. The function fzero uses a combination of bisection, secant, and inverse quadratic 
interpolation methods. 

In Figure 2a we present the steady state solution of equation (10) which reads 

In Figure 2b we present the solution of (10) using the same stress ramp as in the exper- 
iments presented in Ref. [34]. Clearly, for intermediate values of r, the decreasing stress 
branch of the solid concentration lags behind the increasing branch and a hysteresis is 
clearly visible. Having now both r and $ we can calculate 7 from (11) using (5) for nij). 
The values for m, ^^ and e are obtained from the fit performed in the rheometer. We 
obtain the elastic modulus G from a linear fit of the experimentally measured dependence 
r = G7 using the experimental data. The parameters K^^ K^ and w were chosen such 
that the numerical flow curve is in close agreement with the experimental flow curve. 
In Figure 3a we compare the model discussed above with the experimental data. An 
excellent agreement is obtained, particularly at the point where 7 = 0. 

Finally, we test the model against to oscillatory flow measurements where a harmonic 
forcing r = tq sin(27rKt) and the strain response 7 = 7(i) is monitored. This step is 
needed for a proper validation of the model, as the model involves three adjustable 
parameters {Kd, K^ and w). Note that 7 = d'j/dt. As before, having r we solve for 
$ using (10), then 7 is found using (11), and finally we can use standard numerical 
integration to find 7. The results are presented in Figure 3b. 

One can conclude that using the values of the adjustable parameters that have validated 
the flow curves presented in Fig. 3a a good agreement with the oscillatory measurements 
is obtained without any additional parametric adjustments. Fig. 3b. We consider this 
result as a validation of the model initially proposed in [34] and revisited in Sec. 2 which 
further justifies employing it in the linear stability analysis presented through the rest of 
the paper. 

2.3. Non-dimensionalization 

To conclude this section we present the non-dimensional set of equations. Let 

X = xL, U = Wf/max, t = tL/t/max, P = PpU^a.^, T = TflUnrnx/L 

where ^ = ^s+ ^00 is the total viscosity of the completely yielded fluid. Substituting all 
these into (l)-(4) we get: 
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Figure 2: (a) Steady <&. (b) Transient <&. Both with Ty = 6.5Pa K^ = K^. = 0.3 and 
w = d.l. 
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Figure 3: (a) Flow curve, (b) Strain time series. In bot plots we have m = 3, fic 
0.0102, e = 0.01, Ty = 6.5, Kd = K.r = 0.3 and w = 0.1. 



11 



| + (i.V)« = -Vp+|LV.4 + JjV.f (13) 

Vu = 0, (14) 

V 

f + VFe/i (7) <I> f = fi.{j)j, (15) 

5 + («-V)<l> = -5(f)a> + />)(!-<!>). (16) 
at 

where jUj. is the ratio of the solvent viscosity to the total viscosity, i.e. ^s//i. Finally: 

A(7) = I (1 - ^r)^/'" + I r:^ I I , (17) 

5(f) = Kjl + tanhr^^jY (18) 

/(f) = ie, A-tanhT^^^jy (19) 
where 

The Reynolds (Re), Weissenberg (We) and Bingham (B) numbers are defined: 



i?e= """ — , M^e= -■--""— , 5 



/i L tyiiiax/^ 

with the relaxation time Xjj = Hoo/G. 

Through the rest of the paper only non-dimensional variables will be used and the 
"hat" shall be dropped for simplicity. 

3. Linear stability analysis for shear-thinning regularized 
Casson's Model 

First we consider the linear stability of the regularized Casson model. We note that our 
main goal is to investigate the effect the sohd and solid-fluid regions have on the stability 
of the flow. It is known that a viscoplastic fluid is stable to all infinitesimal perturbations 
(at least no instabilities have been found numerically) see [27] , therefore a natural second 
step towards our goal is to ask ourselves how the presence of a highly stratified viscosity 
affects the stability of the flow, thus we consider the case <1> = 0. Linear stability analysis 
for a regularized Bingham model was carried out by Frigaard and Nouar in [11]. There 
they showed that as e ^ the velocity field using a regularized viscosity converges to 
the velocity field using the non-regularized viscoplastic viscosity. While considering the 
problem of linear stability they showed that the spectrum for the regularized problem had 
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some physical spurious eigenvalues, not present in the Bingham problem. The stability 
of the eigenvalues depended on e, we will discuss further about this in Sec. 3.4.3. Here 
we are not interested in the distinguished limit e ^^ 0. Based on experimental data we 
are aware of that typically O(10~^) < e < 0(10^^), therefore for the rest of the paper 
we fix e = 0.002 unless otherwise stated. 

For the rest of this section we set <I> = 0, /Ug = and fi = ^^, then the equation of 
motion reads 

^._V,+ J-V.r, (20) 

along with the incompressibility condition V • w = 0. 
The constitutive equation reduces to: 

r = /i(7)7, 

with ;u(7) satisfying (17) with /i^ = 0. 

3.1. Baseflow profile 

We consider shear flows of the form U = {U{y), 0, 0) with y € [—1, 1], then (13) reduces 
to: 

^ = Re^ (21) 

ay ax 

with the usual non-slip boundary conditions U{ — 1) = U{1) = at the walls. It is clear 
that in this situation 7 = \dU/dy\ and because U{y) reaches its maximum at the center 
of the channel we have: 7 = —dU/dy if y € [0, 1] and 7 = dU/dy otherwise. Having this 
in mind and integrating once (21) we have the following algebraic equation for 7: 

Iy|^i?e + M7)7 = 0. (22) 

dx 

We solve this non-linear equation for 7 using fzero in MATLAB for each point of our 
discrete domain. Once we have the approximation for 7 we integrate it with respect to 
y and use our boundary conditions to get U{y). We feed this result to an user-built 
function in MATLAB and use fzero again to find Re{dp/dx) such that maxj^({7(y)) = 1. 

We present different examples of the base flow for different Bingham numbers in Figure 
4, clearly as we increase B we can see the existence of a pseudo-plug region around the 
center line of the channel. Recall that this is not a true plug, i.e. the velocity there is 
not constant. 

3.2. Linearised momentum equation and tangent viscosity 

As is common in linear stability analysis we consider an infinitesimal perturbation (eu', ep') 
superimposed upon the base flow and linearize the momentum equation (20) around 
{U,p) to get: 
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Figure 4: Base flow proflle for the regularized shear-thinning Casson model with e 
0.002, m = 2.5 and B = 0, 1, 2, 3, 4, 5. 



du' 

— + {u'- V)U + {U ■ V)u' = -Vp + V • r', 

where r' is the stress perturbation given by: 

where jl is the viscosity perturbation and it is given by: 



(23) 



(24) 



(25) 



Given the fact that the flow we consider is unidirectional, it can be shown that (see 
[25]): 



where 



_ KU)jij{u') for ij^xy,yx 
''^ ^ fJ't{U)jij{u') for ij=xy,yx, 

d-yxy 



(26) 



(27) 



is the tangent viscosity. For a one-dimensional shear flow, the tangent viscosity is deflned 
by lit = d.Txy/d'jxy whereas the effective viscosity is defined as /x = Txy/jxy In Figure 5 
we show some examples of the non-linear viscosity n(y) and /ii(y) for increasing Bingham 
numbers, note that the functions are smooth, n > fit and that maxy /u(y) = maxy /ii(y) ~ 
B/e. For a detailed description of the tangent viscosity concept we refer the reader to [25]. 
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Figure 5: (a) Non-linear viscosity /i(y) and (b) tangent viscosity fitiv) with e = 0.002, 
m = 2.5 and B = 0, 1, 2, 3, 4, 5. 
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3.3. Modified Orr-Sommerfeld and Squire equations 

Writing equation (23) in terms of the normal velocity v'{x,t) and the normal vorticity 
ri{x,t) = du'/dz — dw' /dx and assuming modal solutions of the form: 

v'{x, t) = v{y) exp[i(ax + I3z — ojt)], 
r]{x, t) = fj{y) exp[i{ax + (3z — ujt)], 

we get the following eigenvalue problem for the frequency lo: 

ri)(:;)^Kr)' 

where 

Re 

+ ^{D' + k')[{f,t - f^KD' + k% (29) 

Ci = -^{D' + k')[if,t-^)D], (30) 



Re A;2 

' Re k^' 



C2 = ^{DU)-^D[{f,t-fi){D' + k% (31) 



S = aU + ^^fiA + ^{Dij)D + ^^^D[{fit-fi)D], (32) 

with k'^ = a^ + 13'^ , D = d/dy, A = D^ — k'^. Together with the boundary conditions 

v = Dv = r] = at y = ±1. (33) 

3.4. Bounds for the Squire and Orr-Sommerfeld modes and 
one-dimensional stability 

In this section we present bounds for the eigenvalues of the modified Orr-Sommerfeld and 
the modified Squire operators, equations (29) and (32). We follow closely the work by 
Joseph [18] and Davies and Reid [8]. For the rest of the section instead of working with 
the frequency oj as the eigenvalue we consider the wave speed, c = oj/a unless otherwise 
stated. 

As in [18] we define 

r2_ /" u,{n)|2. 



and take v (z H, where ^ is a complex- valued Hilbert space completed under the norm 
/f by the addition of limit points of sequences of functions in C^([— 1, 1]) satisfying (33). 
The following isoperimetric inequalities hold in Ti: 
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Lemma 3.1. Let v ^ H, then 



2 

4 






X2 
^^2 



v^^ 



X2 
^^3 



(2.365)^ 
where A3 is i/ie smallest eigenvalue of a vibrating rod with displacement v satisfying (33). 



3.4.1. Bounds for the Squire and Orr-Sommerfeld modes 

Theorem 3.2. (Damped Squire modes) Let c{a, /3,Re) be any eigenvalue of the homo- 
geneous modified Squire's equation 



1 



—Srj = c rj, 
a 



Then 



r?(±l) 



0. 



Ci < -P' 



aRe 



(34) 



where p = min{^t(y)[ — 1 < y < !}■ 

See A.l for the proof. 
Theorem 3.3. Let c{a, Re) be any eigenvalue of the modified Orr-Sommerfeld equation 



a 



-Cv = ciD"^ - a'^)v 



v{±l) = v'{±l) = 0. 



(35) 



Let q = ioaax{\U'{y)\ : —1 < y < 1} and p = nim{fitiy)\ — 1 < 2/ < 1}, then we have the 
following results: 



1. 



a < 



7r2(7r2+a2) 



-\-a' 



2a aRe \ vr^ + Aa"^ 
2. No amplified disturbances (modes with Ci > 0) of (35) exist if 

qaRe < p ■ f{a) = max[Mi, M2], 
Ml = Agvr + 2^/2^3 



(36) 



(37) 



M2 = AsTT + ttq; , 



where A3 = (2.356)^ 



3. 



U" ■ < • 



uL„ <o<u, 



u„ 



< Cr <U, 



max ■ '-'mtn 



u„ 



+ 



2U"- 



U" < 



'-'min ~r 



7r2 + 4a2 
2U" ■ 

^ mm 
7r2 + 4^2 



+ 

max + 



2U" 

^ max 



2U" 



7r2 + 4a2 



(38) 



< Cr <U, 



max- 
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See A. 2 for a proof. 

We should note that these results are not restrictive to regularized viscoplastic fluids, 
they hold for any shear-thinning (we are using the fact that ^u is a decreasing function of 
7) viscosity function. 

3.4.2. One-dimensional stability 

From Theorem 3.2 we have that the Squire modes are always damped, thus for the rest 
of this section we consider only the modified Orr-Sommerfeld equation and a = 0. If 
instead of working with the perturbation of the normal velocity v we consider a stream 
function formulation, i.e. u = ^y and v = —^x, and taking a modal solution of the form 
^ = (j){y) exp[i(ax — wt)], then we can interchange v for (j) in (28). Note that (p G Ti and 
Theorem 3.3 holds for this formulation. Therefore using (89) we have a bound for the 
imaginary part of the frequency oj: 

qahlp _ {ReyMll + 2aHl + aHl)) 
'"''^Il + a^ll n + a^P, ■ ^^^' 

Here we are interested in one-dimensional perturbations, i.e. a = 0. This is a special 
case in terms of boundary conditions. The fact that v = does not mean that = 
at y = ±1. This boundary condition can be obtained if the perturbation to the flow 
rate j[_-^ u{y)dy is taken to be zero and a uniform pressure gradient in the x-direction is 
allowed, see [36]. Now, even for the case a = we have (j) ^71 and (39) becomes: 

a;. < -^, (40) 

where we used the isoperimetric inequalities defined in Lemma 3.1. Therefore, one- 
dimensional perturbations are always linearly stable. Note again that this result is valid 
for any decreasing function ^(7). 

3.4.3. Solution of modified Orr-Sommerfeld equation and code validation 

There is no equivalent of Squire's theorem for generalized Newtonian fluids but some 
authors have performed several numerical tests for large range of axial and transverse 
wave numbers (a and (3 respectively) for different non-linear viscosity functions ^(7), 
see [24-27]. Their numerical results show that the lowest critical Reynolds number is 
obtained for spanwise homogeneous perturbations (/? = 0, a 7^ 0). We should also note 
that in the situation of a homogeneous streamwise perturbation (a = 0, /3 7^ 0), the 
imaginary part of u satisfies: 

1 j\p^{^p^\Dv\'' + \D^v + PM')<^V , 
'^i = — -;^ i < 0) (41) 

Re J^j2?«|2+^2|^|2dy 

which comes from letting a = 0, multiplying by the complex conjugate of v, say v* , 
integrating by parts and taking imaginary part in (28). Note that if a = then C 
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and S decouple, thus the flow is unconditionally stable. Using these facts we would 
only consider 2D perturbations in the streamwise direction for the rest of the section. 
We know from Theorem 3.2 that all the Squire modes are always damped, thus we are 
concerned only with the modified Orr-Sommerfeld equation: 

ac{D^ - a'^)v =a[U{D^ - a^) - D^U]v - ^^D{fiD)v 

Re (42) 

Together with the boundary conditions 

V = Dv = at y = ±1. (43) 

We solve the eigenvalue problem (42)- (43) by discretizing using Chebyshev polynomials 
in the usual fashion, as described for example in [39]. In order to validate our code, we 
benchmark by solving the Newtonian problem m = 1 and B = 0, i.e. n = fit = ^ 
in (42) and compare with results obtained by Mack [21]. In Figure 6a we present the 
first 33 modes calculated using our code with N = 150 and the ones reported by Mack 
with the eigenvalues of the two computations overlaying one another. Figure 6b shows 
the error norm between our calculations and Mack's for iV = 50 and N = 150. As 
expected the least stable eigenvalue has converged but the split of the S*- branch is present 
when N is not large enough. In order to stably compute the spectrum of (42)- (43) we 
need substantially more nodes than for the Newtonian problem. In Figure 7 we present 
the calculations using N = 150, 200, 250. The first thing to note is that even with 
A^ = 150 the S-branch of the spectrum has not converged. Another important fact about 
the spectrum is the presence of an R-branch, which was first documented by Frigaard 
and Nouar in [11]. They found that this branch is physically spurious, meaning that 
these eigenvalues do not exist for the Orr-Sommerfeld operator corresponding to the 
viscoplastic model and their stability depends on e. As e ^ we can find unstable 
modes from this R-branch even for moderate Reynolds numbers. We have found that 
by choosing a smooth regularized non-linear viscosity (either the one due to Becovier 
et. al. or the one due to Papanastasiu) increases the range of e for which the R-branch 
remains stable and the least stable eigenvalue belongs to the A-branch as expected. We 
have fixed e = 0.002 for all our results and we have checked that the critical eigenvalue 
is indeed in the A-branch. For the rest of the section we fix the number of nodes to be 
A^ = 250. 

3.5. Numerical results 

Stability results for the regularized Casson model are presented in the following section. 
First we introduce the "plasticity number" defined as: 

PI = BRe = ^. 
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Figure 6: (a) Spectrum for Newtonian case, present code with N = 150 (o), Mack [21] 
(+). (b) Error norm of present code with A^ = 50 (diamonds) and A^ = 150 
(squares) vs Mack [21]. AU for Re = 10000 and a = 1. 
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Figure 7: Spectrum for regularized shear-thinning Casson's model, N = 150 (o), A^ = 200 
(+), and N = 250(0). For Re = 10000 and a = 1. Ah calculations with 
e = 0.002. 

Note that PI only depends on the rheological properties of the fluid and the geometry of 
the problem, thus as we increase Re in our analysis PI will remained fixed. All the results 
that follow show the critical Reynolds number for the regularized model normalized with 
respect to Re Newt = 5772.2. In Figure 8a we fix tti = 2.5 and with the increase PI we 
enhance the stability. These results are not surprising, we are mimicking a viscoplastic 
fluid which is known, at least numerically, to be linearly stable [27] with a pseudo- 
plastic fluid with a very large zero-shear rate viscosity. The use of this type of non-linear 
viscosity has also been shown to enhance stability [25] which is also represented in Figure 
9a where we fix PI = 10000 and vary the power-law index m. Figures 8b and 9b show 
the corresponding critical wave number for increasing plasticity number and increasing 
power law index respectively. Note that these values remain in close proximity to a = 1.02 
which is the critical wave number for the Newtonian case. 

4. Linear stability analysis for Elastoviscoplastic model 

In this section we study the linear stability for the elastoviscoplatic model presented in 
equations (13)- (16). 

4.1. Baseflow profile 

For the base flow we consider unidirectional shear flow of the form U = {U{y), 0, 0), thus 
equations (13)- (16) reduce to: 
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Figure 8: (a) Normalized critical Re for increasing PI and m = 2.5. (b) Critical wave 
number . All calculations with e = 0.002. 
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Figure 9: (a) Normalized critical Re for varying power law index m with PI = 10000. 
(b) Critical wave number. All calculations with e = 0.002. 
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dU\ , dT^y ^ dp 



dy \ ^ dy J dy dx 

T^^ = 2Wefii\U'\)<^U'T^y (45) 

T,y = fii\U'\)U' (46) 

Tyy = (47) 

$ = 1^ (48) 

where T = Jl/2 {T^.^ + 2T^y). Substituting (46) into (44) we get: 

i.(„, + ,,ic/'i))f)=ii.| m 

This equation is solved in exactly the same way as (21) in Section 3.1. Once we have 
7 = \U'\ from (49) we solve the following nonlinear equation for Txx'- 

{f{T) + g{T))Txx - 2We {i^{\U'\)uf /(T) = (50) 

where we have used the definition of Txy in (46). We now solve equation (50) using fzero 
in MATLAB® for each point y in our discrete domain. Finally we can construct <l> using 
(48). 

In Figure 10 we present the solution of equations (44)- (48) for increasing Bingham 
numbers and fixed Weissenberg number We = 1. As we can see the velocity profile 
presents the pseudo plastic behaviour which in this case, in addition of a intensely strati- 
fied viscosity, we have a viscoelastic plug, due to the presence of non-zero normal stresses 
inside this region. We should note the existence of a thin region {0{w)) where both solid 
(viscoelastic) and fluid phases coexist as it's clearly seen in the plot of $«. 

4.2. Orr-Sommerfeld equation for Elastoviscoplastic model 

As before, we consider an infinitesimal perturbation {eu' , ep' , er, e^') superimposed onto 
the base flow and the field equations are linearized around (U, P, T, <l>), thus we have: 

^ + (w' • V)U + (U ■ V)u' = -Vp' + ^V • 7(it') + -^V • r' (51) 

ot Re Re 

V • w' = 0, (52) 

r' + Wei^fiCf' + t') + ($/i + $V) T ) = m(«'), (53) 



'dt 
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Figure 10: Base flow profile for the elastoviscoplastic model for increasing Bingham num- 
ber B. From top left to bottom right: velocity, U{y); shear stress Txy] normal 
stress, Txx and steady state structural variable $5. Parameters for this calcu- 
lation: e = 0.002, We = l,m = 2.5, Kr = K^ = 0.3, w = 0.1 and ^^ = 0.1. 
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As there is no version of Squire Theorem for this type of models for simphcity we restrict 
ourselves to two-dimensional perturbations below. We introduce a stream function ^, 
such that: u' = ^y, v' = —^x- We consider modal linear disturbances of the form: 

^ = 0(y)e'("--*), r4. = r,,(2/)e'(""-"*), ij = xx, xy, yx, yy, cl>' = ^(y)eH"— *) , . 

(55) 
Denoting D = -r-, the eigenvalue problem for lo is: 



dy 

1 



Re 



+-^[iaD{Txx - Tyy) + (D^ + a^)Txy] (56) 

Ke 

[ujWe^^iTxx = [1 + mWe^^iU] Txx - 2We^ixDUTxy - iaWe^nDTxx(t> - 2We^fj,TxyD'^(j) 
-2ia [We<^fiTxx + n]D(l)- 2<^nDUTxy^ (57) 

iujWe^fXTxy = [1 + iaWe^fiU] Txy - We^nDllTyy - fitD'^cp 

- [iaWe^fiDTxy + We^fia'^Txx + a^^t] (p (58) 

iojWe^fiTyy = [1 + iaWe^fiU] Tyy + 2\afiD<f> - 2We^fia^Txy<p. (59) 

iioip = [aD<^cl)+[[aU + g{T) + f{T)]ip + [{Kr-Ka)^-Kr]F{y)sign{DU)Txy 
+ [{Kr - Kd)$' - K,.$] Fiy)Wefi\DU\Txx, (60) 



where 



l_tanh2^^ , , 

F(y) = ^ "' J= 61 

w^2(WefiDU<l>)^ + 1 



and with boundary conditions: 

(l) = D(t) = 0, at y = ±1. (62) 

As it is customarily done in linear stability analysis for viscoelastic flows we introduce 
the elasticity number defined as: 

Re' 
which is the ratio of kinematic viscosity and relaxational diffusivity: ei depends only on 
the properties of the fiuid and the flow geometry. For the rest of the paper we replace 
We by elRe. 

4.2.1. One-dimensional stability 

In this section we consider one-dimensional disturbances of (56)-(60). In order for (j) 
to satisfy boundary conditions (62) we proceed in the same way as for the regularized 
viscoplastic fluid, by letting the perturbation to the flow rate J_^ u{y)dy to be zero and 
we impose a uniform pressure gradient in the x-direction. Letting a = in (56)- (60) we 
are left with the following eigenvalue problem: 
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D'^Tc.y = -iReujD'^cl)- UrD'^cp, (63) 

(1 - ieeUe^fiUj) t^^ = 2etRe^ ixDU T^y + 2elRe^ iiT^yD'^ <p + 2^fiDUT^yip, (64) 

{l-[eiRe<^fiu)T,,y = ^xtD'^ct). (65) 

{l-'\etRe^lxuj)Tyy = 0, (66) 

(ia; - g{T) - f{T)) ^ = [{Kr - Kd)^ - Kr] F{y)s\gii{DU)T^y 

+ [{Kr - Kd)^^ - Kr^] F{y)emefi\DU\T^^ (67) 

We should point out that from equations (64)- (67) a continuous spectrum exists 
and this consists of purely imaginary eigenvalues in the strip [—l/eiReinax{fi), —oo) 
(note that for our choice of regularization parameter e, we always have (K^ + Kr) > 
l/€iRem.ax{n)). One has to be careful when considering this part of the spectrum, note 
that as <& — 7> then Wj -^ — oo , clearly these eigenvalues are always damped, thus we 
consider only the case Tyy = 0. Also note that (64) and(67) decouple from (63) and (65), 
then we multiply (63) by the conjugate of (j), say </>*, (65) by r* and integrate by parts 
to get: 

/ Ta;yD^(j)*dy = [Reu / {Dcpl^dy -fir \D^(t>\^dy, (68) 



/ {I - e£Re^iJ.uj)\Txy\'^dy = fitD'^(l)T*ydy. 



(69) 



Expanding in real and imaginary parts and applying the mean value theorem for 
integrals as necessary we have the following equation for the imaginary part of oj, 

uji = i i < (70) 

fit,Rej_^ \D4>\My + elRe^li J_^ \T^y\^dy 

where 

^ti=^t(6), ^ = ^(6), and 11 = ^{^^) for some 6,6,6 ^ [-1, !]• 
Therefore, (56)- (60) is linearly stable to one-dimensional perturbations. 

4.2.2. Code validation 

As a first test for our code we consider the case jir = 0.5, $ = /U = /Uj = 1 and 99 = 
in (56)- (60), by doing this we recover the Oldroyd-B model and compare our results 
with the work by Sureshkumar and Beris in [40]. We set Re = 3960, We = 3.96 (or 
el = 10~^), a = 1.15 and have used N = 128 and N = 256 to solve for the spectrum. We 
show the results in Figure 11, the least stable eigenvalue is a; = 0.3409 -|- 1.9888 x 10^'^i, 
the value reported in [40] is a; = 0.3409 -|- 1.9696 x 10~^i. We have convergence of the 
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discrete spectrum and the difference of the least stable eigenvalue between computations 
is 2.7 X 10-12^ 

We now turn our attention to the full model (56)- (60). For the rest of the section 
we consider the wave speed c = w/q as our eigenvalue. The first thing to note is the 
location of the continuous spectrum, from (57)- (60). For simplicity let us define the 
Deborah function as De{y) = ei^^. We can see that the real part of the continuous 
spectrum is in the line Cr = U and the imaginary is in the line q G [a, oo), where: 

a = max{2i^d, 2Kr, WeB/e}. 

Clearly, the continuous spectrum is always stable but one has to be careful when doing 
computations. Note that for y in the pseudo-plug region (where T < B), De{y) is very 
large (recall that fi ~ B/e), once we begin to approach the yield limit De{y) decays 
to zero very fast as <l> — > and ^ ^ 1. This means that the continuous spectrum will 
approach zero when y is on the pseudo-plug region and it will be pulled down very fast 
towards minus infinity as T approaches B. This could cause severe numerical problems, 
$ is actually zero when the fiuid is fully yielded. We discretize the problem using the 
same technique as in Section 3.4.3. In order to prove the convergence of our code we 
compare the spectrum of a regularized version of (56)- (60) with the one of the full model 
using two different eigenvalue solvers in MATLAB®. Let 

*''^^^~\ ^ for <P{y)<^, ^ > 

where i? <C 1 is a small parameter. We replace <& with $r in De{y) and solve the general- 
ized eigenvalue problem using the function eig in MATLAB®, we choose i? such that we 
get finite results. We then solve the full problem (56)- (60) using eigs in MATLAB®, this 
function provides the reverse communication required by the Fortran library ARPACK 
which is based upon an algorithmic variant of the Arnoldi process called the Implicitly 
Restarted Arnoldi Method. This function can calculate k eigenvalues based on a user 
defined a, thus we choose k such that we get finite results and a is taken to be the 
least stable eigenvalue of the regularized Casson model calculated in Section 3.4.3. We 
present the results of these calculations with N = 250 in Figure 12a, the spectra of both 
problems overlap each other. Note the presence of the almost undisturbed spectrum of 
the generalized Newtonian problem (^-branch, P-branch, S*- branch and ii- branch), the 
continuous spectrum corresponding to equation (60) centered at Cj = {Kd + Kr)/a'^ (we 
have assumed K^ = K^ for this example) . Because the choice of parameters (see caption 
in figure), the pseudo-plug region is very small we have only a few eigenvalues corre- 
sponding to Cr = U "^ 1 and q = —l/a^De{y)Re. Note that there exists a set of discrete 
eigenvalues which emanates from the region close to —\/a'^De{^)Re which distorts the 
A, P, R branches of the modified Orr-Sommerfeld operator. The difference on the least 
stable eigenvalue between the problem using ^^ and the one using <I> is 5.1101 x 10^^. 
This result is not surprising due to the fact that this eigenvalue belongs to the wall 
modes (A-branch), there the fluid is fully yielded and it does not feel the effect of elas- 
ticity, when we use <I>r close to the wall we can assume that e£ '^ 0(10^^) for this choice 



28 



-0.2 



-0.4 



-0.6 




o N=128 
o N=256 



0.2 0.4 0.6 0.8 1 1.2 



Figure 11: Code validation. Spectrum for Oldroyd-B fluid, N = 128 (o), N = 256 (D). 
For Re = 3960, e£ = 10"^ {We = 3.96), /3 = 0.5 and a = 1.15. 

of parameters. In view of these results we use the function eigs for the full discretized 
problem (56)-(60). In Figure 12b we present convergence results for A^ = 250, 350, 450, 
the change on the least stable eigenvalue is on the sixth significant figure. For the rest 
of the paper we fix A^ = 450 unless otherwise stated. 

4.3. Numerical results 

Next we present the results of our cacluations. Figure 13a shows the critical Reynolds 
number for the elastoviscoplastic model as a function of plasticity number, PL As with 
the regularized Casson model the existence of a pseudo-plug region (stratified viscosity) 
is suffciente to greatly enhance the stability of the flow. The critical Reynolds number 
appears to be a monotone increasing function of plasticity number, just as with the 
regularized viscoplastic model. Clearly, in Figure 13b we can see that the behaviour of 
the critical wave number Oc is very similar to the wave number of the Casson model. 
We should also note that the inclusion of a highly viscous viscoelastic fluid as a plug 
destabilises the flow when comparing it with the regularized model (shown as the dotted 
curve in Figure 13a). In relative terms, when PI = 1000 the critical Reynolds number 
for the elastoviscoplastic model is 2.66% smaller than the critical Reynolds number for 
regularized model. For PI = 10^ this percentage increases to around 6%, this is due to 
the fact that the pseudo-plug and solid-fluid regions increase and are closer to the wall. 

As a way to understand better the effects of the elastoviscoplasticity on the instability 
mechanism we will turn our attention to the analysis of the distribution of the production 
and dissipation of disturbance kinetic energy across the channel. The linearized Reynolds- 
Orr equation is derived by multiplying equation (56) (equation (42) for the regularized 
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Figure 12: Code validation for elastoviscoplatic model, a) Spectrum for full model (o) 
and for regularized model using ^r (+) with A^ = 250. b) Spectrum for full 
model with N = 250 (o), N = 350 {+), and N = 450(0). For Re = 10000 
and Q = 1. Parameters for these calculations: e = 0.002, e£ = 10~^, m = 2.5, 
Kr = Kd = 0.3, w = 0.1 and fi^ = 0.1. 
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Figure 13: (a) Normalized critical Re for increasing PL (b) Critical wave number. Pa- 
rameters for these calculations: e = 0.002, ei = 10~^, m = 2.5, K^ = K^ = 
0.3, w = 0.1 and /U^ = 0.1. 
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model) by the conjugate of the stream function perturbation (or velocity perturbation 
in (42)) and integrating over a periodic domain. This can be written as: 

%^^fe>-ife) (T2, 

where j =regularized (reg) or elastoviscoplastic (ev) and (•) = J_i-dy . The left hand 
side of (72) is the temporal variation of the averaged kinetic energy, the first term in the 
right hand side is the exchange of energy between the base flow and the perturbation 
and the last term in the equation is the rate of energy dissipation. Explicitly we have for 
the regularized model 

h,re9 = \Dv\'^ + a^\vf (73) 

l2,reg = a[DUregiVrDVi - ViDVr)] (74) 

/3^^g^ = 4a^fi\Dv\'^ + fit\D'^v + a^v\'^ (75) 

and for the elastoviscoplastic model 

h^ev = |Dvp + a>P (76) 

l2,ev = a[DUev{(prDct)i - (PiDcPr)] (77) 

h,ev = fir\D'^4' + (y.'^4>? 

+ «[('^ra ~ Txx)rD(l)i - {Tyy - Txx)iD4>r\ 

+ T^yAD^(t}r+a'^^r)+T^y,{D^(l)i+a'^^i)- (78) 

Following Govindarajan et al. [13], we compare the normalized space-averaged energy 
production, 

and normalized space-averaged energy dissipation 



Reihj) 



TJ = ^Jf , . (80) 



At criticality these two quantities balance. Positive values of I2J indicate where in the 
flow domain energy is supplied from the base flow to the perturbed flow. It should be 
mentioned that even though that the energy production and energy dissipation balance 
at critical conditions, the mechanism of instability is governed by l2^i and not by /sj-. 
We refer the reader to [13] for an extensive discussion about Reynolds stress distribution 
analysis. In Figure 14a we can see how F^ and T~ balance at critical Reynolds number 
Rcreg = 13364.41 for PI = 1000. The inclusion of an elastic plug an its effect on the 
stability of the flow is clearly seen in Figure 14b. We keep the same Reynolds number as 
for the regularized case. Note that both F^ and F^ are in the same orders of magnitude 
as before, but the balance is slightly broken. The energy production has increased but 
note that the energy dissipation has gone negative close to the centre line of the channel 
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(insert figure) whereas for the regularized model is always positive due to the presence 
of a highly viscous pseudo-plug. This is due to the presence of elastic forces. Even 
though the elastic plug is confined to a very small region away from the wall, it is enough 
to break the balance between energy production and dissipation and stability is lost. 
Experimental evidence of the existence of an elastic plug can be found in [28] . 

Now we explore the effects that increasing elasticity will have on the stability of the 
flow. In Figure 15a we present the normalized critical Reynolds number for increasing ei. 
Note that the change is minimal, the same happens for the critical wave number shown 
in Figure 15b. This is not surprising either, elastic forces are confined to regions of very 
low shear-rates where the value of the nonlinear viscosity is very large, thus an increase 
on ei does not make any difference. This is also seen in Figure 17a, where e£ = 10~^ and 
r+ and rg;^ are almost unchanged with respect to the case ei = 10^^. In Figure 16a-b 
we show the effects of increasing the parameter w which in turns makes the solid-fluid 
region wider. Note that for small values of w things remain more less unchanged but as 
soon as w gets closer to 1 we see a steep decrease on the critical Reynolds number. This 
again is not too surprising, by increasing w the region where solid-fluid coexist stretches 
towards the wall. When it; = 1 we have that both <I> and Txx are different from zero at 
the wall. This is clearly seen in Figure 17b where the normalized production of energy 
r+ has increased in the critical layer. 

5. Discussion 

We have studied linear stability analysis of plane Poiseuille flow of an elastoviscoplastic 
fluid. Our results show that there is an increase in the critical Reynolds number for 
two-dimensional perturbations as the yield stress increases. The main difference between 
our study and previous works is that we consider a rather new yielding scenario recently 
suggested in [34] according to which the transition from solid to fluid regime is not direct 
but mediated by a solid-fluid phase coexistence regime within which the material behaves 
as a viscoelastic fluid. Thus, we were interested in understanding how the presence of 
elasticity modifies the hydrodynamic stability of the flow. The relevance of a particular 
yielding scenario to the stability problem might not be obvious at a first glance as the 
yielding typically occurs at low Reynolds numbers (the experiments described in [34] were 
performed at Re<l) whereas the loss of hydrodynamic stability emerges at significantly 
higher Re. However, it has been demonstrated experimentally in [14] that the transition 
to turbulence in the pipe fiow of a viscoplastic fluid (Carbopol 940) is simultaneous with 
the breakdown of the unyielded plug initially located around the centerline of the pipe. 
These results suggested us that the loss of the hydrodynamic stability and the yielding 
are in fact connected, which motivated us to revisit the stability problem in the context 
of the yielding scenario suggested in [34] . These results are in a close agreement with the 
ones of a regularized viscoplastic model, this is not surprising due to the fact that the 
viscoelastic core is confined to a region away from the wall. We have also shown that the 
existence of a region where solid and fiuid structure coexist destabilizes the flow. 

We have chosen the regularized Casson model for the viscosity function but we suspect 
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Figure 14: Energy budget for Re = 13364.41 and PI = 1000. (a) Regularized Casson 
model, r+ = r~ = 7.152 x 10"^ (b) Elastoviscoplastic model, r+ = 7.351 x 
10^^, r^ = 7.033 X 10^'^. Parameters for these calculations: e = 0.002, 
ei = IQ-^, m = 2.5, Kr = Kd = 0.3, w = 0.1 and ^r = 0.1. 
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Figure 15: (a) Normalized critical Re with PI = 1000 and increasing ei. (b) Critical 
wave number. Parameters for these calculations: e = 0.002, m = 2.5, Kr = 
Kd = 0.3, w = 0.1 and fi^ = 0.1. 
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Figure 16: (a) Normalized critical Re with PI = 1000 and increasing w. (b) Critical wave 
number. Parameters for these calculations: e = 0.002, e£ = 10~^, m = 2.5, 
Kr = Kd = 0.3 and /Xr = 0.1. 
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Figure 17: Energy budget for Re = 13364.41 and PI = 1000. (a)Elastoviscoplastic model 
with el = 10^3, r+ = 7.356 x 10"^ F" = 7.03 x W^ (b) Elastoviscoplastic 
model with tt; = 1, r+ = 7.436 x 10"^, F" = 7.037 x 10"^. Parameters for 
these calculations: e = 0.002, m = 2.5, Kr = K^ = 0.3 and fir = 0.1. 
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that the stabihty behaviour is quahtatively similar for other models such as the regu- 
larized Hershel-Bulkley and Bingham models, and therefore that our results are quite 
generic. 

There are practical limitations to be acknowledged. First, the domain is a geometric 
idealisation of a planar geometry of large aspect ratio. Depending on the actual geometry 
it becomes necessary to consider other effects, e.g. entry/start-up effects of the flow, cur- 
vature or imperfections of the walls, etc. On the other hand in the experimental setting, 
the point at which instability is actually observed is very sensitive to control of apparatus 
imperfections and the level of flow perturbations. For example, in Hagen-Poiseuille flow of 
Newtonian fluids one typically observes transition to turbulence starting for Re > 2000. 
However, an experimental flow loop in Manchester UK produces stable laminar flows for 
Re ^ 24, 000, [17,29], and stable flows have even been reported up to Re « 100, 000, [30]. 
This suggests that enhanced stability may be achieved experimentally, where predicted 
by the linear theory. 

A. Bounds for Orr-Somerfeld and Squire modes 

A.l. Proof of Theorem 3.2 



Proof. The homogeneous modified Squire equation is obtained by setting t; = in (28). 
This gives 



CT] =Ur] + -!—ij.{D'^ - k^)ri + -^DfiDr] 
aRe aRe 



(81) 



aRe k"^ 

Multiplying (81) by the conjugate rj* , integrating by parts over y € [—1,1] and using 
boundary conditions we get: 



r?|2 dy= / U\7]\'^ dy + 



aRe 



M|I?r?|2-fc2H2) + ^(^,-^)|D^|2 



dy. (82) 



Using the second mean value theorem for integrals and taking real and imaginary parts, 
we get: 

cril = U{yi)ll (83) 



Cil 



-1 



J-'O 



aRe 



(84) 



with yi G [—1,1] to be the mean values. The first bound immediately comes from the 
fact that Umin < U{yi) < Umax- Notc that for shear thinning fluids, /i(yi) > min^u > 
min^j > 0. We also have iJ-tiVi) > min^i. Thus, for p = min^t > 0, 



Q/q < 



1 



aRe 
_P_ 
aRe 



[il + eil] < 0, 



(85) 
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so all the Squire modes are damped since Jq is positive. Dividing through by /q and 
applying our isoperimetric inequality Ii/Iq > 7r^/4 given in Lemma 3.1 we get our result 
as required. D 

A. 2. Proof of Theorem 3.3 

Proof. Using the fact that the Squire modes are always damped, we multiply (42) by the 
conjugate v*, integrate by parts, use boundary conditions, take real and imaginary parts 
and use the second mean value theorem for integrals to get: 

Ci{lf + a^/^) = {Q- Q*) - {aReyHa'^illf - {aRe)-^fit f \v" + a^vp dy, (86) 

Cr{lf + a2/2) = f [U\v'\^ + (a^U + U"/2)\v\^] dy, (87) 

where Q = {i/2)fU'v*'v dy, fi = ^(yi), fit = l^t{y2) with yi, y2 G [-1,1]- We can 
expand the last integral of the right hand side of (86) and integrate by parts. This gives 

Ci{lf + a2/2) = (Q - Q*) - {aRe)-\ia^fLlf - jit{ll - 2a'^lf - a^I^)]. (88) 

We will now proceed with the proof. 

As before, note that jl > min^u > min^Ui, thus we divide by {if + a^ Iq) , using Cauchy- 
Schwartz inequality on the first term of the right hand side of (88) and our definitions 
for p and q, we have 



qhlo {aRe)-\^a^pll + p(/| - 2a^ll - a^/^)) 



II + a2/2 II + a2 j2 

qhlp {aRe)-^p{ll + 2a^/f + g^/g)) 

II + a2/2 /2 + a2/2 



(89) 



Clearly, equations (87) and (89) are the same as equations (3) and (4a, b) in [18] respec- 
tively, with the inclusion of the positive constant p in (89). Therefore the proof follows 
in the same way as in [18] and our results hold. D 
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